function [likel,ssvec,flag_ok,lht] = likelCumReg(parest,parvec,parposest, funcmod, Y, trainvec, solveopt, addsol)
% =========================================================================
% likelCumReg.m 
% Formerly known as Kfilter_tagtv_general.m 
% Computes likelihood of model with Cumulator but no Split
%
% This is a file for the LMJ project 
% It computes the likelihood a DGSE model where there are time aggregeation
% issues and the model is written in terms of cumulator variables.  
% ========================================================================
%% 
%% DESCRIPTION 
%
%% I. State space representation ot time invariant model 
%
% $$ y_{t} = Z_s s_{t}+ C_sZ $$
% 
% $$ s_{t+1} = T_{t} s_{t} + R_{t} \eta_{t}$$
% 
% $$ V(\eta_{t}) = Q =SDX'SDX $$
%
%% II. Inputs 
% 
% *PAREST*: vector of parameters 
% 
% *FUNCMOD*: model function handle 
% 
% *Y*: [nobs nz] data vector, with NANS 
% 
% *TRAINVEC:* [2 1] vector indicating which row to start in and finish for
% the evaluation of the likelihood. Note: if empty [1 nobs] will be used. 
% 
% *SOLVEOPT*: setting for solution, e.g. for NL solver 
% 
% *ADDSOL*: additional input needed to solve model 
% 
%% III. Construction of time aggregated model 
%
% $$ y_{t} = Z_{t}\alpha_{t} + CZ $$
%
% $$ \alpha_{t} = T_{t}\alpha_{t-1} + R_{t}\eta_{t} $$
% 
% where $$ \alpha_{t}=[s_{t};M_{t}] $$ 
% and $$ M_{t} $$ is a cummulator whose transition equation depends on whether
% aggregating SUMS or AVERAGES 
% 
% $$ M_{t}=\Psi_{t} M_{t} + As_{t} $$
%  
% Since cumulator placed at the bottom of the state, the model solution
% must supply the additional output FIELD, a structure supplied by the solution of 
% the model,that contains structures used to set-up extended state-space 
%
% .NZM      Number of m series 
%
% .NZQ      Number of q series 
%
% .A        [NZQ NS] matrix for the transition of the cumulator 
%            Selects the rows of the state with the elements of the
%           cumulator 
%
% Treatment of the constant changed in July 21 2010. Now the extra rows of
% C for the expanded state space will be created by A*C, hence enter in C
% the constant of those observations that will be part of the cumulator.
% Field CADD no longer needed. 
%
% .FLAG_TAGSUM      == 1 if aggregation is for sums 
%                   == 0 if aggregation is for averages 
% 
%% IV. Construction of time varying matrices 
% Current code defines 3 type of matrices 
% $$ G_{t}, R_{t} $$ that vary deterministacally with 
% month j=1,2,3. 
% 
% The observation matrix $$ Z_{t}=W_{t}*Z $$, where $$ W_{t} $$ is
% extracted every period by removing NANs from the data. 
% Note, however, that the code assumes observations are either observed
% every period or every 3 periods. 
% 
% 
%% V. Output 
% 
% *stt*: [nobs ns] vector of one-sided states 
%
% *etamat*: [nobs nx] vector of smooth disturbances 
%
% *smooth_st*: [nobs ns] vector of smooth states 
% 
% *yfor*: [nobs nz] vector of one step ahead forecasts for observables. 
% Note: notation is not exactly as in data unless NANs all grouped at the 
% end of the observation vector
% 
% *vfor*: [nobs nz] vector of one step ahead forecast errors for observables. 
% Note: same as yfor, notation is not exactly as in data unless NANs all grouped at the 
% end of the observation vector
%
% *logL*: [1 1] log likelihood with constant included 
% 
% *countmat*: [nobs ns nx] decomposition of all states when 1 shock at a
% time turned on. 1 shock per page. 
%
% *counobs*: [nobs nz nx] decomposition of all obserbvables when 1 shock at a
% time turned on. 1 shock per page. 
% 
% *countszero*: [ns nx] decomposition of state S(1) when 1 shock at a time
% turned on. 1 shock per column. 
% 
% *likelvec*: [nobs 1] vector of likelihoods 
%  
%% VI. Additional Notes 
% a) It is assumed the number of periods is divisible by 3 and that the first period 
% is month 1, such that the quarterly series are not observed. 
%
% b) Data are demeaned inside the code  
% 
%% VII. Related codes 
%
% *KFILTER_TAGIN.m* Assumes a deterministic pattern in the 
% variation of Z, from M to Q. Requires additional outputs from the DSGE. 
%
% *KFILTER_TAGIN_FLEX.m* Does not require deterministic pattern in
% variation of Z. Works with any model with NANS and otherwise time
% invariant matrices. 
%
% File Created on 06/25/2010
% Modified        07/01/2010
% 7/21/2010      Changed treatment of constant 
% ========================================================================
%% 
%% ELEMENTS OF THE CODE 
%% Initialize parameters
likel=-1e40; flag_ok=0;
parvec(parposest)=parest;
%% 1. Solution of the model
%     FIELD contains structures used to set-up extended state-space 
[G,R,C,eu,SDX,Z,structOne,ssvec,structTwo]=feval(funcmod,parvec,solveopt,addsol);
if ~isequal(eu,[1;1]);return;end

% Note: Above change momentarily for debugging of kfilterCumSplit: RAB
% should remove junk,junk from output if there...


% =========================================================================
%% 2. Set-up extended states space 
%     See buildState.m by AB for an alternative 

%% 2.1 Dimensions and time invariant matrices 
[T, ny]=size(Y); 
ns       =size(G,1); 

if size(Z,1)~=structOne.nzHF;error('Rows in Z must match number of variables observed in High frequency'); end 
nx  = size(SDX,2);

Ztilda  = [Z zeros(structOne.nzHF,structOne.nzLF); ...
    zeros(structOne.nzLF,ns) eye(structOne.nzLF)];
% Selection matrix 
A   = structOne.A;
Ctilda   = [C; structOne.A*C]; 
% Number of periods of High frequency per 1 in Low frequency 

 % type of temporal aggregation
if structOne.flag_tagsum == 1
    sigma =[0 ones(1,structOne.NUMIT-1)];
else
    sigma =(1:structOne.NUMIT);
end

%% 2.2 Page matrices with pages corresponding to quarters {1,2,3,4} 
Tpag  = zeros(ns+structOne.nzLF,ns+structOne.nzLF,structOne.NUMIT);
Rpag  = zeros(ns+structOne.nzLF, nx ,structOne.NUMIT);
for ii=1:structOne.NUMIT 
    if structOne.flag_tagsum == 1
        Tpag(:,:,ii) = [G zeros(ns,structOne.nzLF); ...
            structOne.A*G eye(structOne.nzLF)*sigma(ii)];
        Rpag(:,:,ii) = [R; A*R];
    else
        Tpag(:,:,ii) = [G zeros(ns,structOne.nzLF); ...
            (1/sigma(ii))*structOne.A*G eye(structOne.nzLF)*(sigma(ii)-1)/sigma(ii)];
        Rpag(:,:,ii) = [R; (1/sigma(ii))*structOne.A*R];
    end 
end 

% =========================================================================
%% 3. Demean data      
if any(C~=0)==1;Y=Y-repmat((Ztilda*Ctilda)',[T 1]);end;
Y=Y';

% =========================================================================
%% 4. Initialization
%  Begin with Tpag(:,:,1), Rpag(:,:,1) which should be ZERO for the cummulator 
shat=zeros(size(Ztilda,2),1);
pshat=lyapunov_symm(Tpag(:,:,1),Rpag(:,:,1)*(SDX'*SDX)*(Rpag(:,:,1)'));
    
%% 5. Forward filter 

% 5.1 Storage matrices 
lht      = zeros(T,1);
casevec  = kron( ones(T/structOne.NUMIT,1) , (1:structOne.NUMIT)' );
Zdim    = zeros(T,1); 
W        = eye(ny);

for ii=1:T;
   
    % 5.2 Remove NANS from ytt
    % Construct $$ Z_{t}=W_{t}Z $$
    ytt=Y(:,ii);
    
     % Determine W and position of the NAN  
    ind =~isnan(ytt);
    
    % create Z*,y* from
    Ztt = W((ind==1),:)*Z;
    ytt = ytt(ind);
    Zdim(ii)=length(ytt);
    
    % 5.3 Kalman routine
    [shat,pshat,lht(ii)]=kf_dk(ytt,Ztt,shat,pshat,Tpag(:,:, casevec(ii) ),...
        Rpag(:,:,casevec(ii))*SDX' );
end


%% Likelihood Calculation

lht     = lht(trainvec(1):trainvec(2));
likel   = -0.5*(sum(Zdim))*log(2*pi)+sum(lht); 
flag_ok = 1;

%% End of File
end